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ABSTRACT 


The problem of potential subsonic flow around 
complex configurations is considered. This problem requires 
the solution of an integral equation relating the values of 
the potential on the surface of the body to the values of the 
normal derivative, which is known from the boundary conditions. 
The surface of the body is divided into small (hyperboloidal 
quadrilateral) surface elements, , which are described in 
terms of the Cartesian components of the four corner points. 

The values of the potential (and its normal derivative) within 
each element is assumed to be constant and equal to its value 
at the centroid of the element. This yields a set of linear 
algebraic equations. The coefficients of the equation are given 
by source and doublet integrals over the surface elements , . ■ 

Closed form evaluations of the integrals are presented. The 
results obtained with the above formulation are compared with 
existing analytical and experimental results. 
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SECTION I 


FORMULATION OF THE PROBLEM 
1, 1 Introduction 

A general theory for compressible unsteady potential aero- 
dynamic flow around lifting bodies having arbitrary shapes 
and motions is given in Refs. 1 and 2. Application to finite- 
thickness steady and oscillating wings in subsonic flow is 
given in Refs. 3, 4, and 5. A general numerical formulation 
for complex configuration? In subsonic flow is considered in 
Ref. 6. Applications of this formulation are considered here. 
For simplicity, only the incompressible flow is considered 
in-detail*. In this case, the problem is governed by the 
Laplace equation with prescribed normal derivative on the 
body [exterior Neumann problem for the Laplace equation) with 
an additional complication due to the presence of the wake 
Cof unknown geometry) . 

The problem of the evaluation of the steady, incompressible 
potential aerodynamic flow around an aircraft of arbitrary 
configuration can be analyzed by solving the integral equation 

tsv v - ^^3°^ <1-1 

_ . S' 25 

where o is a surface surrounding the aircraft and the wake ' 

(Fig. 1) . For the moment, it will be assumed that the wake 

does not exist. The effect of the wake is considered in 

Section III. 


* Subsonic oscillatory flow is considered in Appendix A, 


The value of is obtained from the boundary condition 

(tangency condition) 


(x.fV)~U a (n^^L^o 

T>K on ' ^on J 


Cl. 2) 


or 


rop 


- ~n x =-yt'£ 


^ n (1.3) 

The integral equation can be studied by dividing the surface 
CT into N small finite elements G'c (see Fig. 2) to yield 


11 ”-* ~+fr,.v(±)]d^ 

6* 

Applying the mean value theorem one obtains 


Cl. 4) 


(1.5) 


**<$ ' " u **< % 

where % is a suitable mean value of inside the element 

which will be approximated by the value of ^ at the centroid 


M 


P ** of the element, (T^ . 

By satisfying Eq, (1.5) at the centroid, P v ‘ , of the 
element 6^ , (h - 1, 2...N) yields 


jff Yjfn.v(jr)d^ (kzl,l,),...N) (i.e> 
lti h~t & k fh * 


3 


where r^ is the distance of the centroid of the element 6 h 
from the dummy point of integration in the element 
Equation Cl. 6} is equivalent to* 


where 


and 


with 


C ^h-k - f %} = ( K } 
c /^ s i? If 

* 



(1.7) 


( 1 . 8 ) 


Cl. 9) 


a. io) 


1.2 Surface Geometry 

Let the geometry of the element <5** he described by 


p {~f[ -S’) 


( 1 . 11 ) 


where ^ 1 and -j? v are the generalized curvilinear coordinate. 
Then the two base vectors Cu are given by {Fig. 3) 

-> *5 p 

The unit normal to the surface is given by 


( 1 . 12 ) 


— > 
n 




(1.13) 


a, * a. 


The effect of the wake is not considered here ' (see Section 
III) . 


4 


and is directed according to the right-hand rule (Fig. 3.). 
The surface element c}<5" is given by (Fig. 3 ) 




d-f ' x a* d-f ’ * | 5, x | cJ-f ! d“f s 


(1.14) 


1.3 Expressions for b^ and c 


'hk 


Combining Egs. (1.10), (1.13) and (1.14) yields 


'h-k = 


2TT 


JJ 


a/ xa 2 l 


i i V^- 


«/ f a ¥ 


(1.15) 


Similarly, combining Egs. (1.8), (1.13) and (1.14) yields 


where 


c 


1 


H = 2tt ).j S; i X n x * v (y)<W ' d-f * 





(1.16) 


(1.17) 


and 


r~ 1 ?.) 


(1.18) 
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SECTION II 

HYPERBOLO IDAL ELEMENT 


2 . 1 Introduction 

Consider the equations 


^ + $ 1 ’'? * + fi -f 7 

* ^ c + Jt"? i- 7 + Jj -J’ 7 . 


( 2 . 1 ) 


or r in vector notations 

P = Pc + ?.-f + R9 - My (2,2) 

This represents a hyperboloid. The lines ^ = const and = 
const are clearly straight lines. Consider the hyperboloidal 
element defined by the above equation with 


-1 < -J < 1 

- 1 < 7 ± ) 

(2.3) 

The centroid of the element is P c ( ^ = 0)* The corner 

points of this element are 

“ pc r Pi + + fj (“?= + ! , y.r? 

P.- = P. : -r f, " R " Pj (-?=<!. V = -i) 

;fc - Pc - R +R-p s (-f *- i' , •)-*') 

£- = R, - P, - P, * p, (-f .-i, 5 =-i ) 


The inverse relation is 


Jl , 

?c -- -r ( P*. f i- t p.,+ P„) 



ijj 

f. = i ( ?« ■ * i - - P-) 



'iH 

3i . 

U.J 

= + 1- i) 



CJ 

i { Ri - P- - it ■ f ?-) 


£2.5) 

Jli 

Note that the four boundaries of the 

element 


!!. £ 
iili 

are straight lines given by 



LiU 


-1 4 ^ 

■s 1 

r ~, — i 

l! 

LUi 

P-- (jp«- p,) + tf-fDj- 

7 

£ 1 

\ 

;|j 

L_J 

h %)$■■ 

-i -i 

i 1 

— i 

i 

j 

P * (?c - p.? v 

- 1 <-f 

t 


( 2 . 6 ) 


Next, assume that the surface of the aircraft is divided 
into curved quadrilateral elements with four corner points 

p ++ , p + _ , p_ + , p . Then, as mentioned in Section X, these 

elements can be replaced by the hyper bo loida 1 element (des- 
cribed above) which goes through the four corner points p . . , 

p + _ t p_ +/ p (see Fig. 4 ) . it may be noted that the surface 

is continuous since adjacent elements have in common the 
•straight: line connecting the two .common corner points . It 


may be noted also the P c is the centroid of the hyperholoidal 
element and hence it will be indicated as 




{2.6) 


2.2 Geometry of Hyperboloid Element 

The geometric quantities introduced in Section I can be 
written for the hyperboloid element described above. Letting 


¥ ~ F j 


Equation (1.12) yields 


Ft? 




JW.I 


'D % 


This yields 


-- a, =( p, * P 3 j/)x(p z x p 3 j) 

Pi *■ Pz t' Pi * P 3 F + B *'Pz f 

since pi^ x p 3 — 0.. 

Note tha,t r with- present notations- 

l = r k ~ p‘-? a, = ? a -P fh +aT+l?+'BF? 


Po + p, f -? a?+p 3 


where 


(2.7) 

(2.8-a) 
(2 . 8-b) 


(2.9) 


C2 t i0L 


is the .vector connecting the centroid P - of the element t 
to the one r P ^ t of the element , 


2.3 Expressions for end 

Introducing the functions I g and I Q (indefinite source and 
doublet integrals) such that 


r 5 _ l | a { * a.z\ 


'dJQy 2 .U r 


(2.12) 


and 








Zlt 


(2.13) 


Equations Cl. 15) and (1.16) may be rewritten as 


b kl = it 0,0- Is f -<)- h (-■/, i ) - r 5 (- 1,-0 


(2.14) 


/J— Ij* f/i ~/J - I]) (-/; /) -+ Xif (~{,-0 


(2.15) 


In Ref, 6 , the expressions for I g and I Q have been obtained as 
(the proof is included here in Appendix b). 


iui.ih *&■>'{ If fcfefe-) 


(2.16) 


r s 


■ ' -<l W|- M '■ 


u 


N 


^ |ll!a<| + i'tti U 


p=-7-& |5||«4+ i-n bm( Iff 'I a 1 - — ) 

- UJ : V r Vkk-a.xa^ / 


g_* cC|xai 

(2.17) 




•■T.r 




2.4 Solution of Problem 

The wake contribution is shown in detail in Section III , 

Equation 1.7 can now be solved for % , (see also Ref. 7 . 

* 

The linearized Bernoulli Theorem will further give the pres- 
sure distribution: 



- ^ Z 1% 

3 x 


(2.18) 


2.5 Wing and Body Geometry 

This subsection presents the geometry of the wings and 
wing-body combinations which were used for obtaining the 
results presented in Section IV. 


2,5,1 Wing Geometry 

First/ the wing geometry is presented. The projection 
of the wing on the xy plane is given by 


" X.e. {$) 

where L.E. and T.E. stand for leading 
respectively. At' any point along the 


(2.19) 

edge and trailing edge, 
span, b, the chord is 


given by 


cfj) r X- t - fa) - o) (2.20) 


The planform described in Eg . (2 .19 ) can be transformed 

into a rectangular one by the transformation. 




K - * 1 , 1 , 
*r. e . “ 

r 


{ 2 . 21 ) 


- 10 - 


Fdr all the cases considered, the wing is symmetric so that 

E = t H [Z . (2.22) 

where h is the thickness of the wing. 

Combining Eqs . (2.21) and (2.22) , we can rewrite the wing 
geometry as ^ 

X = fV-iS-e.) 7+ 



Z = t h/2 

with 


(2.23) 


c ^ f ^ i 
J 


(2.24) 


When the wing has an angle of attack, . cO, the position of any 
point on the surface is changed to (see Fig. 5). 


A = t ? Sin 

J = 5 

£ 5 E Co5 - XSin d 

(2 

For small ©( and small thickness ratio, t? / Eq. (2.25) can be. 
approximated as 


X - X 

0-3 ■ - (2.26) 


" 11 " 


The surface of the wing is divided into small quadrilateral 
elements. The corner coordinates of the wing boxes are stored 
and used as input for generating the hyperbolordal surface 
elements. 

A further transformation of the wing geometry is necessary , 
in order to obtain a more accurate evaluation of T near the 
leading edge and near the tip of the wing, where • t P. varies 
more rapidly than in other areas. The transformation is given 
by 


f -- T* 

and it changes the uniform mesh of size 


(2.27) 


At = C./NS3C 

AO =• U/2WY _ 

(where NX is the number of boxes in the p-direction and NY 


(2.28) 


is the number of boxes in the' 5 ! "direction along the semispan) 
into a nonuniform one, as follows: 

r along the -direction, the boxes are larger near the 
wing root and smaller near the tip; 
along the f -direction, the boxes are smaller near 
the leading edge and increase in size as we approach 
the trailing edge. 

2. 5 , 2 Thlchness distribution 

For most cases considered here the thickness is taken to be 

■ v-k- VC m * £5! T)('- (2.29) 

2 - 
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where C vnct* is the. maximum chord, "ZT is the thickness ratio 

__ fl max 

(2.30) 

where h M <i x is the maximum value of the thickness for a given section 
Another wing profile investigated is the circular biconvex one. 

The planform is described by Eg. (2.26) and the thickness distri- 
bution is given by 


4 = 2 


c +h 




HI a* 


4 h 


MAX J 


-U- 


^LE ~%T£ 


) - 


t\ k 




2 4 (2.31) 

As before, c represents the chord and h m a*- represents the maximum 
thickness. The radius of the circular arc is 

c + h m <\ x _ 

(2.32) 


R = 


4 

Once the wing profile is selected, a variety of wing planforms 
(rectangular, swept, delta) can be easily generated by simply 
changing the leading and trailing edge angles of the wing. 


2.4.3 Body geometry 

Eor the wing-body configuration, the wing is attached to a 
body of revolution composed of a forebody of length Lp and radius 

r -2- 

(2.33) 

a midsection of length and constant radius r = i/Z . and an 

l (Fig. 6) . The 
surface of the body is also divided in boxes. Care must be 
taken in matching the coordinates of the wing and the body boxes at 
the wing root, in order to prevent discontinuities in surface. 


aftbody of length L_^ and constant radius f - 
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SECTION III 
THE WAKE 

3,1 Dynamics of the Wake 

As mentioned in Section J, th.e surface <f in Eq. 1,1, 
surrounds the body and the wake. The effect of the wake, 
disregarded in Section X, yields an additional terr.i in Eq, 

(1.4) , given by 2 ' 5 

Iw=h-JJ (3 - 1) 

with tr 

A*? ~ ^ (3.2) 

This represents a distribution of doublets with intensity 
A^. The geometry of the wake is not known. An iterative 
procedure can be used to solve the problem: consider the 

surface of the wake divided into small elements. Assume 
initially that the wake is composed of straight vortex lines 
(see next subsection) : then find the values of c Pi and then 

evaluate the velocity at the corner of the elements. Find 
a new location for the corner of the element such that the 
elements approximate the stream surface emanating from the 
trailing edge and repeat the procedure mentioned above. 

However, a simplified treatment of the wake is considered 
here, 

3,2 Simplified Treatment of the Wake 

The simplified treatment of the wake used here consists 
of assuming that the wake is composed of straight vortex^ 


-14- 


lines emanating from the trailing edge and parallel to the 
x-axis (direction of the flow)* For this case, the surface 

A 

of the wake is divided into infinitely long elements, , 

with two edges parallel to the x-hxis. These elements are the 
continuation of the elements of the wing having an edge in con- 
tact with the trailing edge (Fig* 7 ) , 

Hence, by assuming that (in view of the KUtta condition 
that ^(&¥)/c)x~0 at the trailing edge), the value of 
can be approximated by the value at the centroid of the element 

A -r* 

% the contribution J-v/ (see Eq. 3.1) is given by 

% (3.3) 

with 

r X* = ± 2^Jj Ku * ^ (3,4) 

fit 

for the elements with an edge in contact with the trailing 
edge, and • 

^k - 0 (3.5) 

for the others. 

In order to evaluate the integral in Eq. (3.4), it is con- 

ft 

Venient to consider that the element Ojj, is the limit of the 

parallelepipedal element obtained by truncating the element 
* 

tfs, at the finite distance (Fig. 8 ) . The limit is obtained 
by letting 

5^5 * y.1 ( 3 . 6 ) 
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go to infinity; note that u = i since two’ edges are parallel 
to the x-axis. Note that (see Fig. 8 ) 


if- j? = = p m 

?i - pi = ■£, 

B . . ft . 5T, 

P 3 - O 


(3. 


(3. 


It is convenient to separate the contribution from the 
trailing edge (T - r t ) and the edge that goes to infinity 


(T ' i): 


iw - + iTrfit - Jw C _ 1 (3 

where (note that ^ + ^i) 

5 s • *-i * Q-i)- Ay^ C Pn-rf • A n ^ ) (3 

with 

f/nj = Pm ■+ Pcf (3 

while (note that ^ (|*= i) - •+ 2^2 ) 

P ; ' 4 “ ^ 111 lv s -i x ^i /^i 

w7 -(l**)-CE*5) \ 


7) 

8 ) 

.9) 

. 10 ) 
. 11 ) 


(3.12) 


and similarly (note that ^ (Jr - t) s 




" fcoAl 


' r (t£x 

\ l* M* • a. *«: 




I 




f (1 **)' (|_X Pet,) 

l'| I If'* 


/ 


T 5 - 1 

T= 


- 1 

[ R>>(j x j 

M 

[ Pmj * 

&) 

1 ^<0 1 

1 i ?u ■ 

X X frf 




~ f f*>«* ’ jj*£ ) ( /^ • Ptf) ~ ( Pjhrf • Pcf )( [m«r X ) 
\ [ PmJ P*<f * A X ft* | 


(3.13) 


/ 


SECTION IV 


RESULTS 

4 . 1 Introduction 

The formulation outlined in the previous section has been 
'translated into a computer program. In this section , the 
numerical results obtained for wings and wing-body combination 
in steady and oscillatory subsonic flow are presented and 
compared with existing ones. 

4 . 2 Wing Results 

Figures g. and 10- show the general distribution of the 
potential <j> and of the lift coefficient C L for a wing in sub- 
sonic flow. 

Figures 1 la and 1 2 present a comparison with the experimental 
and analytical results of Lessing, Troutman and Menees (Ref. 

11). The results are for a wing of aspect ratio AR ■= 3. The 
wing profile is a biconvex circular arc, 5% thick, with sharp 
leading and trailing edges. Figure Ha shows the thickness 
effect (the pressure distribution on the upper and lower surfaces 
of the wing) for the wing at zero angle of attack and Mach 
number M - .24. Figure 12 represents the lift distribution 
f or = 5° and M = .24. The whole wing was divided into 196 

elements, or NX = NY = 7 for the upper right hand side of the 
wing in the X and Y directions, respectively. A convergence 
study for the problem given for Fig. 11a is presented in Fig. 
llb t The curves plotted in Fig. xib are the distributions 
of the velocity potential along / =0 for different numbers 

of elements. It is shown that the results obtained by using 
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100 elements on the whole wing, or NX = NY = 5 are accurate 
enough for general analysis. 

Figure 1 3 shows the lift distribution per unit angle of 
attack, Cjtot , for a rectangular wing with AR = 1 and M = . 2 , 
These results, obtained with NX = NY - 7 are compared to the 
ones of Cunningham and Kulakowski and Haskell (Refs. 12 and 
13} . 

Figure 1 4 shows the distribution of CJ?<* for a tapered 
swept wing with aspect ratio AR =3, taper ratio TR = .5, 

■Ai - 45 ° and M = . 8 . The results obtained with NX = NY 
= 7 are compared with the ones of Cunningham ( Ref. 12) and 
Kolbe and Bolts ( Ref. 14), 

Figure 1 5 shows the distribution of the section lift 
coefficient per unit angle of attack, t for & rectangular 

wing with aspect ratio AR = 4 and M - .507. The results, 
obtained with NX = NY = 7 and then 10, are compared with the 
ones by Yates (Ref. 1ET ) . 

4.3 Wing-Body Results 

Results for wing-body configurations in steady subsonic 

flow are compared in Fig. 1 6 a with the results presented by 

17 

Labrujere, Loeve and Slooff. The results were obtained for 
M = 0 and a rectangular iidpositioned wing with chord c = 1 , 
span b = 6 , thickness ratio X = .9, and - 6 °. The body is 
at zero angle of attack and is composed of a forebody with 
length = 2 and 

r = o.s - o.izF (K - x,,.f (5.1) 
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a midsection of length — 1 and radius, r = .5, and an aft- 

body of constant radius r = . 5 and length - 9, No wake is 
used on the body. The number of elements is 200 on the whole 
configuration C NX = 5, NY - 4 on the wing, NX = 2, NY - 3 on 
the forebody, NX = 5, NY = 3 on the midsection and NX - NY = 

3 on the aftbody) . In Fig. 1 6a, the distribution of the wing 
section lift coefficient is presented. In Fig. lfh the dis- 
tribution of the difference of the velocity potentials on 
the upper and lower surface of the fuselage are shown. These 
results are related to three circumferential stations as shown 
in the figures. Results for subsonic oscillatory flow are 
presented in Figs. 17 and 18 for the same wing considered 
in Fig. 12, oscillating in bending mode 
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with k = cu c/Z t-L, = .47, and M - .24 (NX = NY = 7). Figure 
17 presents a comparison with the results of Lessing, et al., 
while in Fig. 18, the results are compared with the ones by 
Lessing, et al., and Albano and Rodden (Ref. 15). 
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SECTION V 
CONCLUDING HE MARIAS 

A finite -element method for linearised steady and oscilla- 
tory subsonic potential aerodynamics around complex configurations 
has been presented. Numerical results have been obtained for 
wing-bodv configurations in steady flows and for finite thickness 
wings in oscillatory flows. However, the formulation is general 
can be applied to more complex configurations, although, in this 
case, the accuracy of the method remains to be assessed. For 
instance, in the case of wing-tail interaction, the roll-up of 
the wake might play an important role. A rolled-up wake geometry 
.is being considered for lifting surfaces. 

The method is very flexible and simple to use; the use of 
quadrilateral hyperboloidal elements . (which can be used to yield 
any arbitrary closed surface) , defined in terms of their corner 
points, is one of the original features of the method. Another 
original feature of the method is the simplicity of the expressions 
for the coefficients, due to the vector formulation of the problem. 
The method is also accurate and fast , despite the fact that no 
effort has been made yet to minimize the computation time; as an 
example, the computer time for the results for subsonic wing- 
body configurations (Fig. 6 ), with 200 surface elements is 
775 seconds on the IBM 360/ 50 of Boston University's computing 
center. The results obtained are in excellent agreement with 
existing ones. 
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Fig. 7 Treatment of the wake 
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Figurel lb Analysis of Convergence: Potential Distribution, r*. Versus x/c, at y ~ 0, 
for Recfanguiar Wing With Biconvex Section, in Steady Subsonic Flow, for 
AR = 3, t - 0.05, M - 0.24, a = 0°, NX = MY « S, 6, 7. 




Figure 12 The distribution of the lift coefficient/ eg f on a 
symmetric rectangular wing with AR= 3j,T= 0.05 
5°, M - . 24 and NX=NY=7 for comparison 
with results of Ref. ri 
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FtgurelS The distribution of the section lift coefficient per unit 
angle of attack for a rectangular wing with AR = 4, 

M - .507 and NX=NY=7 and 10 for comparison with 
results of Ref. 16 
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Figure 17 The distribution of lift coefficient, , 

for a rectangular wing oscillating in bending mode 
with k =«c/2U„ « .47, M=.24, AR=3,-c= 0.05, 
NW =20, L V/ =2.5c and NX=NY=7 for comparison 
with results of Ref, 11 • 
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Fig. 18. COMPARISON WITH RESULTS OF REFS. H AND 15 FOR 

RECTANGULAR WING OSCILLATING IN BENDING MODE WITH 
Jc =toc/2U, = .47/ M - .24/ AR = 3 and NX = NY = 




APPENDIX A 


SUBSONIC OSCILLATORY FLOW 


A<1 Integral Equation 

In this Appendix, it is shown how the results obtained 
in the main body of this report can be extended to subsonic 
oscillatory flow. Introducing the variables 


X = X z.-x T=fj £L= slL , 1} 

p* t t t P a * 

/v 

and the complex potential (fi such that 

r <a >> . ;n(T+wx)-i 

$(*.#,*): V tlx* f faxz)e I (a . 2) 


the integral equation for the subsonic oscillatory flow is 
given by 
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where Zl surrounds body and wake. 


A* 2 Boundary Condition 

The boundary condition is given by 
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where Cj) and cj) are such that 
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Next, assume that the motion of the surface consists of small 


harmonic oscillations around a rest configuration , that is 
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with 


and 


\zs = 0 ( 1 ) 
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f A. 15) 


it is easy to see that 
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Neglecting the terms which contain € (of order £ ) and 

separating the steady from the oscillatory terms, one obtains 
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Equation (a* 19) reduces to 
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Finally, neglecting terms of order £ in Eg. (ft.. 18) and terms 
3 

of order £ in Eg. (A. 19), one obtains 
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(where the upper [lower] sign holds on the upper [lower] surface), 
one obtains 
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Equation (A* 27) gives the value of dj>/dH to be used in Eg. 
(A. 3) . 


A 3 Pressure Coefficient 


The pressure coefficient is gr^en by the linearized Bernoulli 


theorem as 
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For oscillatory flow, setting 
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APPENDIX B 

vr* <^r^c^r « , v <r tyir «g v. «. ■ 

Xn this Appendix.^ it will be shown by diffentiation that 
the results obtained for the doublet and source integrals are 
valid for any hyperboloidal element, 

'- BrS- Doublet integrals 

Xhe doublet integral is: 

I» = 'Mjl IX£L*Jk) 

ZTt cb.d 


From Eg. {.1*12), it is apparent that 
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Note that 
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This proves Eg. (2.13) 

B.3 1 Source Integral 

Now, we will prove by dif f erentialtion that the result ob- 



tai.ned for the source integral is valid for any quadrilateral 
planar element. In this case, since the normal 
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and using Eq, (B,8). we get 
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Now, combining Eqs. (B,10), (B.13), (B.14), and (B.16), yields 
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